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Abstract 

The many-body diffusion quantum Monte Carlo (DMC) method with twist-averaged boundary conditions 
is used to calculate the ground-state equation of state and the energetics of point defects in fee aluminum, 
using supercells up to 1331 atoms. The DMC equilibrium lattice constant differs from experiment by 0.008 A 
or 0.2%, while the cohesive energy using DMC with backflow wave functions with improved nodal surfaces 
differs by 27 meV. DMC calculated defect formation and migration energies agree with available experimental 
data, except for the nearest-neighbor divacancy, which is found to be energetically unstable in agreement with 
previous density functional theory (DFT) calculations. DMC and DFT calculations of vacancy defects are in 
reasonably close agreement. Self-interstitial formation energies have larger differences between DMC and DFT, 
of up to 0.33eV, at the tetrahedral site. We also computed formation energies of helium interstitial defects where 
energies differed by up to 0.34eV, also at the tetrahedral site. The close agreement with available experiment 
demonstrates that DMC can be used as a predictive method to obtain benchmark energetics of defects in metals. 

PACS numbers: 61.72.J-, 64.30.Ef, 02.70.Ss 
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I. INTRODUCTION 



The mechanical properties of metals are dominated by the formation and migration energies 
of defects. Experimentally it is often difficult to measure desired defect properties directly, so 
it is important to have accurate theoretical approaches to calculate defect properties. The many- 
body diffusion quantum Monte Carlo (DMC) approach 1 is the most accurate method for systems 
with more than w 30 electrons but has not previously been applied to defects in metals. Until 
today, the most successful quantum mechanics-based defect calculations in metals use density 
functional theory (DFT). However, the approximate functionals used (i) lack sufficient specific 
and universal accuracy, (ii) can not be systematically improved, and (iii) there are now many 
approximate functionals to choose from, all giving different results. Thus, DMC results are an 
ideal candidate to wholly replace DFT when benchmark thermodynamic properties are required. 
Furthermore, the relative energies of reference phases (bulk metals and compounds) would be of 
great value in thermodynamic databases and in the subsequent prediction of phase diagrams (e.g. 
CALPHAD and Thermocalc). 

In semiconductors, DMC calculations of the formation energies of point defects and surface 
energies have shown that important corrections to DFT arise when electronic correlations are fully 
taken into account. Deviations in formation energies of more than 1 eV were found in silicorP 
and diamond 3 . Additionally, activation energies of common chemical reactions obtained by DFT 
methods have been shown to differ substantially from benchmark diffusion Monte Carlo (DMC) 
values 4 and quantum chemical results. For metallic systems, the size of the errors in DFT calcu- 
lations of defects is largely unknown, as more accurate benchmark calculations do not currently 
exist. It is highly desirable to demonstrate the feasibility of DMC, with its increased predictive ac- 
curacy, as a replacement of DFT for challenging systems such as metals, particularly as computer 
power increases. 

Aluminum is an ideal starting point for carrying out initial DMC calculations of defects since it 
is one of the simplest metals with a close-packed fee structure that contains no 3d electrons. As a 
result it is considered a prototype material for testing the validity of theoretical calculations. Alu- 
minum is well characterized experimentally so there is an abundance of data available. There have 
been previous quantum Monte Carlo calculations of the bulk properties of aluminum, however, 
these were done using the less accurate variational Monte Carlo 5 6 and the calculations had large 
statistical noise. 



2 



In this paper, we report well converged results for the bulk properties of fee aluminum using 
DMC. We explore a larger range of volumes in order to compare the ground-state equation of 
state calculated with DFT. Our DMC calculations of the defect properties of aluminum include 
the simplest point defect, the vacancy for which numerous DFT 7 9 and experimental resultsP^ 
are available. We compute the nearest-neighbor divacancy binding energy, a defect that DFT 
calculations^^ have found to be unstable. Since this instability is counter to both experimental 
studies^**' and simple bond counting arguments 17 it is important to perform calculations of this 
defect with DMC, a method that unlike DFT, does not rely on approximations for exchange and 
correlation. 

We also examine two other types of defects: First, self-interstitials can arise due to irradiation 
with energetic particles, through plastic deformation, or through their production in thermal equi- 
librium at high temperatures. We have obtained DMC results for the formation energies of the 
(100) -dumbbell, the octahedral, and the tetrahedral self-interstitials. Second, experimental stud- 
ies of irradiated aluminum show the presence of He bubbles^^. The formation and growth of 
helium bubbles can alter a material's mechanical properties through void swelling, embrittlement, 
and surface blistering^ 22 !. Since irradiation, through He implantation or transmutation, gives rise 
to He atoms at substitutional or interstitial lattice sites, the energetics of these types of defects are 
important. Presented here are DMC calculated formation energies for He at the substitutional, oc- 
tahedral and the tetrahedral interstitial sites, which are compared with previous DFT calculations^. 
We demonstrate that twist-averaged boundary conditions offer far superior statistics over DFT cal- 
culated single-particle finite-size corrections for DMC calculations of defects, and also remove a 
dependence on data from other methods. Overall the calculations reported here used several mil- 
lion processor hours, which is affordable on large computer clusters and supercomputers. 

II. BULK ALUMINUM 

The DMC approach^ is a stochastic method for evolving a wave function using the imaginary- 
time Schrodinger equation. In principle, and in contrast to most electronic structure methods, the 
systematic errors can be measured and systematically reduced^ 5 ! In practice, the most significant 
sources of error are (i) the fixed-node and fixed-phase approximations, a variational solution to 
the Fermion sign problem, (ii) adequately sampling the Brillouin zone, which in contrast with 
insulators is a significant problem in metals, and (iii) pseudopotential error and corresponding 
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locality error (although these may be avoided through all-electron calculations^). 

For our DMC calculations we used the CASINO codeP^I, with guiding wave functions formed 
by a product of Slater determinants for up and down spin electrons and a Jastrow correlation 
function. The single-particle orbitals in the determinants were obtained from DFT calculations 
using the generalized gradient approximation (GGA) for the exchange-correlation term since it 
should perform better than the local density approximation (LDA) where the electron distribution 

"71 

shows large spatial variations as it does at a vacancy in aluminum . (For bulk aluminum GGA 
is more accurate than LDA, see Fig. [3}) For GGA we used the Perdew-Burke-Ernzerhof (PBE) 
form 28 rather than the Perdew-Wang-91 (PW91) form 2 ^ since it gives better values 8 for the vacancy 
formation energy. The DFT calculations were performed using the plane-wave PWSCF codeP 
with Troullier-Martins non-local pseudopotentials. For the non-local pseudopotentials we used 
the locality approximation^ in the DMC calculations. The orbitals were evaluated in DMC via 
real-space cubic splines^. For the DMC simulations we used 5280 walkers with a time step of 
0.01 au, which our test calculations revealed gave time-step errors in Table [I] of less than 0.01 eV. 

Calculations of solids using DMC require finite simulation supercells with periodic boundary 
conditions imposed on the Hamiltonian. This leads to finite-size effects that can be divided be- 
tween single-particle and many -body contributions. In a metal the DMC kinetic energy contains 
large single-particle contributions due to the sharp Fermi surface, which impacts our calculations 
in two ways. 

First, in a real metal the number of orbitals with energies below the Fermi level is usually not 
equal to the number of electrons required in the simulation cell. In DFT calculations one can 
use partial occupations of these orbitals to create a closed shell configuration guaranteeing that 
the charge density has the correct symmetry. In a DMC calculation with a guiding wave function 
containing a single determinant for the spin up and for the spin down electrons the use of partial 
occupations is not an option. Gaudoin, et. alP found differences in aluminum as large as 0.1 
eV/atom in variational Monte Carlo total energies depending on the occupations of the orbitals 
at the Fermi level. As our calculations show in Fig. [3] for aluminum an error on the order of 0.1 
eV/atom is too large to accurately discern the equilibrium lattice constant. 

A second complication for DMC calculations of metals arises as the system size is varied, 
which can produce band crossings as energy levels pass through the Fermi level. This can cre- 
ate discontinuous changes in the nodal surface of the guiding wave function as the symmetry of 
the occupied orbitals change, leading to discontinuities in the DMC energy as shown in Fig. [TJ 
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FIG. 1: (Color online) The black circles with error bars are the calculated DMC ground-state energies of 
bulk fee aluminum obtained using twist-averaged boundary conditions with 10x10x10 twists. The red dia- 
monds with eiTor bars are the calculated DMC ground-state energies of bulk fee aluminum obtained without 
using twist-averaged boundary conditions. The DFT corrections E^ MC = E® MC + E^ FT — E® FT were 
added to the results that did not use twist-averaged boundary conditions. Both sets of calculations were 
done using 5x5x5 supercells containing N = 125 atoms. 

For simulation cells containing 64 or 125 atoms we found that single-particle DFT corrections 
E^ MC = E® MC + E^ FT — E® FT were ineffective in removing the large discontinuities we found 
in calculated DMC energies as the volume was varied. However, the use of well-converged twist- 
averaged boundary conditions'^ was effective in producing smooth energy versus volume curves 
at these systems sizes as shown in Fig. [T] Our calculations with 64 and 125 atoms used T-point 
centered grids with 13x13x13 and 10x10x10 twists, respectively. 

With twist-averaged boundary conditions the remaining finite-size effects arise from many- 
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FIG. 2: (Color online) The circles are calculated DMC ground-state energies of fee aluminum with lattice 
constant a=4.046 A using twist-averaged boundary conditions for simulation cells containing 27, 64, 125, 
216 atoms, and 1331 atoms, with a T -point centered grid of 17x17x17, 13x13x13, 10x10x10, 9x9x9, and 
lxlxl twists, respectively. The statistical error bars are smaller than the circles. The dashed line is a guide 
to the eye. The solid red line is a fit to the data using Eq. ([I]) where N, the number of atoms, is 64, 125, and 
216. This fit has a correlation coefficient of -0.9997. 

body contributions 35 that scale as 

E 00 = E N + c/N. (1) 

As shown in Fig.[2]this equation agrees well with our DMC energies using twist-averaged bound- 
ary conditions for simulation cells containing between 27 and 1331 atoms. 

We used Eq. ([T]) to obtain infinite-size extrapolated DMC total energies from calculations using 
64 and 125 atom simulation cells with twist- averaged boundary conditions using a range of lattice 
constants shown in Fig. [3} The energy points were fit to a quartic and a Murnaghan equation of 
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FIG. 3: (Color online) Upper: Comparison of equilibrium lattice constants and cohesive energies for fee 
aluminum calculated using DMC with backflow (4.030(1) A, 3.403(1) eV), using DFT with LDA (3.960 A, 
4.21 eVj^and GGA (PBE) (4.046 A, 3.52 eV), and experiment (4.022 A, 3.43 eV) with zero-point energy 
and finite-temperature effects removecP. Lower: DMC energies for various lattice constants, with a quartic 
and a Murnae harP fit. DMC data was obtained using twist-averaged boundary conditions and extrapolated 
to infinite-sized supercells with single determinant fixed-nodes. A shift of 0.063 eV was applied based on 
DMC calculations using backflow wave functions (see text). 

state 36 . Both fits yielded an equilibrium lattice constant of 4.030(1) A. This compares well with 
the experimental value, 4.022 A, with zero-point energy and finite-temperature effects removecP. 
In contrast DFT calculations with the local density approximation (LDA) and GGA (PBE) yield 
3.960 A and 4.046 A, respectively. 
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For the DMC cohesive energy we initially obtained 3.341(1) eV using the fixed-nodes defined 
by the GGA orbitals, compared with the experimental value of 3.43 eV with zero-point energy and 
finite-temperature effects removed, 4.21 eV using LDA, and 3.52 eV using GGA (PBE). Although 
the fixed-node DMC results with GGA nodes are already the most accurate, to assess the pos- 
sible nodal error we also performed DMC calculations with twist-averaged boundary conditions 
and extrapolation to an infinite-sized supercell using optimized backflow wave functions (with 
backflow transformations that contained electron-electron, electron-nuclei, and electron-electron- 
nuclei terms 37 ) at the optimum lattice constant and also for an isolated atom. Backflow wave 
functions can be substantially more accurate than single determinant non-backflow wave func- 
tions, typically yielding an additional few percent DMC correlation energy in atomic calculations 
and nearly 100 percent of the correlation energy in the homogeneous electron gas^. However, 
since backflow is too expensive to apply routinely, we have used the backflow result, similar to 
how corrections for all-electrons have been applied 26 , to shift our single determinant fixed-node 
energies in Fig. [3] The backflow cohesive energy is 3.403(1) eV. A complete backflow evaluation 
of the lattice constant might further reduce the residual differences from experiment. 

We expect the backflow correction in metallic systems with atomic number higher than alu- 
minum to be at least as large as our computed correction in aluminum. This indicates that to obtain 
cohesive energies to better than O.leV - other than by fortuitous error cancellation - backflow or 
other nodal optimization must be considered. 

We performed additional DMC total energy calculations of bulk aluminum for atomic volumes 
smaller than those shown in Fig. [3] by following the same procedure of using Eq. (1) to obtain 
infinite-size extrapolated DMC total energies from calculations using 64 and 125 atom simulation 
cells with converged twist-averaged boundary conditions. A Murnaghan fit of the energy- volume 
DMC data was used to obtain the pressure for a range of atomic volumes. These calculated DMC 
pressures are shown in the upper part of Fig. [4] along with pressures calculated with DFT using 
GGA (PBE). At this scale the differences between the DMC and GGA pressures are not visible. 
The solid black curve in the lower part of Fig. [4] shows the difference in pressures between GGA 
and DMC. A common procedure for constructing an equation of state of a material at low temper- 
atures is to shift the computed equilibrium lattice constant so that it coincides with the experimen- 
tal equilibrium lattice constant^!. Following a similar procedure of shifting the GGA calculated 
energy-volume curve so that the GGA equilibrium lattice constant of 4.046 A coincides with the 
DMC equilibrium lattice constant of 4.030(1) A yields a pressure difference between the GGA 
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FIG. 4: (Color online) Upper: The red squares show the calculated pressures of fee aluminum using DFT 
with GGA (PBE). The black circles are pressures obtained from a Murnighan fit of total energy DMC calcu- 
lations performed at a range of atomic volumes. At this scale the GGA and DMC data are indistinguishable. 
Lower: The black solid line is the difference in pressure between the GGA and the DMC calculations at the 
same range of atomic volumes. The dashed red line shows the difference in pressures between the GGA and 
the DMC calculations after the GGA energy-volume curve as been shifted so that the equilibrium volumes 
of the GGA and DMC calculations coincide. 

and the DMC equation of states that increases markedly at smaller atomic volumes as shown in 
the red dashed line in Figj4j This demonstrates that the applying a rigid shift to DFT equation of 
state calculations can result in larger errors than if a shift is not applied. 
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III. DEFECTS IN ALUMINUM 



The results of our calculations of point defects are presented in Table [I] The atomic positions 
were taken from complete structure and volume relaxed DFT calculations using GGA at zero 
pressure. The calculated GGA and DMC vacancy formation and migration energy agree with 
experiment. However, GGA no longer agrees with experiment when "surface" correction s^ 39 * 40 * of 
0. 15 eV and 0.05 eV are added to the GGA (PBE) functional producing 0.82 eV and 0.65 eV for the 
vacancy formation and migration energy, respectively. Previous GGA (PBE) calculations 8 of the 
vacancy formation energy without the "surface" correction using 4x4x4 supercells yielded values 
of 0.61 eV and 0.63 eV with norm-conserving and projector augmented-wave pseudopotentials, 
respectively. This difference with our GGA results is likely partially due to finite-size effects since 
we obtained 0.64 eV for the vacancy formation energy using a 4x4x4 supercell. The GGA results 
in Table [[[correspond to calculations using finite-size converged 7x7x7 supercells. Convergence of 
the GGA defect structures was established by computing the energies using 4x4x4, 5x5x5, 6x6x6, 
and 7x7x7 supercells. The DMC results in Table [I] were done using 5x5x5 supercells with V- 
point centered grids of 10x10x10 twists. The finite-size errors for our DMC calculations are likely 
to be small since the largest GGA energy difference among all of the defects comparing 5x5x5 
and 7x7x7 supercells was 0.02 eV. Shown in Table |Il] are GGA and DMC data demonstrating the 
convergence with supercell size for all the defects considered. 

For the nearest-neighbor divacancy with DMC we obtain a negative binding energy, -0.10 eV, 
which implies that two isolated vacancies are energetically preferred to a nearest-neighbor diva- 
cancy. This agrees with previous DFT calculations which also found a similar negative binding 
energy, -0.05 eV using^LDA and -0.08 eV using 21 GGA (PW91). Thus the disagreement be- 
tween previous calculations and the original interpretation of experimental dataU^, which gave 
positive binding, is likely not a result of DFT approximations. Our results are consistent with the 
reinterpretation 7 of the data. 

Experimentally 42 the (100) -dumbbell was found to be the lowest energy self-interstitial in alu- 
minum. Of the self-interstitials investigated we found that the (100) -dumbbell has the lowest 
formation energy. The calculated formation energy was 2.94 eV using DMC and 2.70 eV using 
GGA. Our DMC value agrees with the experimental estimates of 3.CP 1 and 3.2(5) eV^. For the 
(lOO)-dumbbell we obtain a relaxation volume of 2.25, which agrees closely with experimental 
estimates of 1.9(4) and 1.7(4) 42 . Our GGA (PBE) result is larger than a previous GGA (PW91) 
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TABLE I: Energies (enthalpies at zero pressure) for point defects in fee aluminum (eV) calculated using 
GGA (PBE) and DMC. Formation and migration energies for a single vacancy are denoted HI and H™, 



respectively. Shown in braces are the GGA values with "surface" corrections 7 39 added. Formation and 
binding energies for a nearest-neighbor divacancy are denoted h£ v and H\ v , respectively. Hi, Hi and 
h( are the formation energies for the (lOO)-dumbbell, the octahedral, and the tetrahedral self-interstitials, 
respectively. Formation energies for a He impurity at a substitutional, octahedral, and tetrahedral sites are 
denoted Hg,y H^,y and H^ He , t y respectively. Statistical errors bars for all DMC energies are 0.01 
eV. GGA and DMC calculations were done using 7x7x7 and 5x5x5 supercells, respectively, with atomic 
positions taken from complete structure and volume relaxed GGA calculations. 





GGA 


DMC 


Experiment 


H f 


0.67 {0.82} 
0.60 {0.65} 


0.67 
0.60 


0.67(3) 10 , 0.672H, 0.66(2^ 
0.62P, 0.61(3jSI 0.65(6^ 


H f 

h\ v 

h{ 


1.37 
-0.03 
2.70 


1.44 
-0.10 
2.94 


1.17(7)° 
0.17(5 j3D, 0.2(P 
3.CPI 3.2(5jSI 


Hi 


2.91 


3.13 




Hi 


3.23 


3.56 




H f 

H f 
H f 


1.63 
3.26 
3.33 


1.72 
3.58 
3.67 




"Computed using 


; Hi from Ref. QD) and H$ v 


from Ref. ATI 





resulP3 of 2.43 eV. For the self-interstitials we see differences between the calculated DMC and 
GGA formation energies as large as 0.24 eV. 

The DMC formation energy for a He substitutional defect is 1 .72 eV, while the energies for He 
interstitials are larger than 3 eV. The ordering of these energies is consistent with sites with larger 
free volumes having lower energies. Previous DFT calculations 23 using GGA (PW91) obtained 
1.53, 3.18, and 3.20 eV for He at the substitutional, octahedral, and tetrahedral site, respectively. 
Comparing our GGA and DMC calculations we see differences between 0.09 and 0.34 eV for the 
He impurity. Similar to the self-interstitials, the GGA exchange-correlation errors are larger for 
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TABLE II: Energies (enthalpies at zero pressure) for point defects in fee aluminum (eV) calculated us- 
ing density functional theory (DFT) with GGA (PBE) functional and diffusion Monte Carlo (DMC). The 
symbols for the formation and migration energies are the same as those used in Table [TJ Except where 
indicated the statistical errors bars for the DMC energies are 0.01 eV. The GGA and DMC calculations 
were done with atomic positions taken from complete structure and volume relaxed GGA calculations. 
For the GGA calculations we used 4x4x4, 5x5x5, 6x6x6, and 7x7x7 supercells with T-centered grids of 
13x13x13, 10x10x10, 9x9x9, and 8x8x8 kpoints, respectively. For the DMC calculations we used 4x4x4 
and 5x5x5 supercells. For DMC without twist averaging we applied the single-particle DFT corrections 
E^ MC = E% MC + E^ FT - E^ FT . For the DMC calculations with twist averaging we did not apply DFT 
corrections and used T-centered grids with 13x13x13 twists and 10x10x10 twists for the 4x4x4 and 5x5x5 
supercells, respectively. Entries in the table that were not determined are denoted N.D. 



GGA 


DMC with DFT corrections 


DMC with twist averaging 


supercell 


4x4x4 


5x5x5 


6x6x6 


7x7x7 


4x4x4 


5x5x5 


4x4x4 


5x5x5 


Hi 


0.64 





65 


0.67 


0.67 


0.58 


0.81 


0.66 


0.67 


H™ 


0.59 





59 


0.60 


0.60 


0.50 


0.40(3) 


0.54 


0.60 


H f 


1.36 


1 


35 


1.36 


1.37 


1.31 


1.29(3) 


1.47 


1.44 


H 2v 


-0.05 


-0 


05 


-0.02 


-0.03 


-0.15(2) 


0.33(5) 


-0.15 


-0.10 


H f d 


2.82 


2 


72 


2.72 


2.70 


N.D. 


N.D. 


N.D. 


2.94 


H[ 


2.93 


2 


90 


2.94 


2.91 


2.90 


2.75 


3.18 


3.13 


H{ 


3.53 


3 


25 


3.27 


3.23 


3.27 


3.27 


3.60 


3.56 


H f 

n He(s) 
H f 

n He(o) 
H f 


1.61 
3.26 
3.44 


1 

3 
3 


61 
24 
35 


1.62 
3.27 
3.35 


1.63 
3.26 
3.33 


1.41 
3.29 
2.89 


1.98 
3.45 
3.55 


1.64 
3.53 
3.71 


1.72 
3.58 
3.67 



these defects than for the vacancy. 

Relying on DFT corrections to minimize the single-particle finite-size errors instead of twist 
averaging, yielded poorer convergence for defect energies. For example, the difference in vacancy 
formation energies between 4x4x4 and 5x5x5 supercells was 0.23 eV compared with 0.01 eV with 
twist averaging, while the result for divacancy binding was reversed as shown in Table |TTJ. DMC 
calculations in larger 6x6x6 supercells would be an order of magnitude more expensive, and may 
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still be inferior to the twist- averaged results. The use of twist averaging is essential in metals for 
defects and excitations in which the fractional change in the total energy due to the presence of the 
defect or excitation is inversely proportional to the number of atoms in the supercell, i.e., "1/N" 
effects. 

IV. CONCLUSIONS 

In summary, DMC with twist-averaged boundary conditions can be used to obtain an accurate 
equation of state of aluminum. Our DMC results confirm previous DFT calculations that the 
nearest-neighbor divacancy is unstable in aluminum. Our calculated formation and migration 
energies of point defects show excellent agreement with available experiment, demonstrating that 
DMC can be used to obtain benchmark energetics of defects in metals and used as a baseline where 
no experiment is available. 
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